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ly^ I The Hubbard Hamiltonian is investigated by means of a new variational 

trial wave function. The trial wave function includes either intrasite and 
Y ' nearest - neighbor correlations in an explicit form. To calculate density ma- 

c/3 , trices the method of Kikuchi's pseudoensemble is used. The case of half-filled 
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fermionic band is carefully investigated in the limit of a large number of lat- 



i-^ ' tice sites. The ground state energy and correlation functions are determined 

O I for Bethe lattices with z = 2,4 and 6 nearest neighbors. 
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Short-range correlations have been shown to have a considerable influence on properties 
of strongly interacting fermions. In the most reflned way, the problem is stated within the 
framework of the one-band Hubbard model [0-0]. If fermions of spin 1/2 hop to the adjacent 
sites of a lattice only, the Hubbard Hamiltonian has the following form 

where a]^ (o-ja) is the creation (annihilation) operator of a fermion of spin a on the i-th 
lattice site, {ij) denotes a pair of adjacent sites, rii^ = aj^aja-- 

There are known few exact solutions of the Hamiltonian (|l]): for a homogeneous chain 
(ID) 0], so-called Gutzwiller approach (GA) for a lattice of inflnitive dimensions D= oo 
PJ^ and some other special cases @J31- Lattices of intermediate dimensions, especially 
2D and 3D ones, are of great practical importance, and a number of works is devoted to 
their investigation (see |^,^ and reference therein). In particular, the Gutzwiller's trial wave 
function which is exact in the limit of D= oo was applied |TD|-|T^. This wave function entered 
in a numerical procedure of the variational Monte Carlo (VMC) method |TD|| or was used in 



diagrammatic expansions of GA in terms of 1/D [Tl||T3]. In the last case, one had to base 
on a strong assumption that the ground state of 2D and 3D lattices only slightly differed 
from that of the D= cxd lattice. Actually, the problem of lattices of flnite dimensions lies 
in necessity of inclusion of spacial correlations, especially short-range ones. As was shown 



phenomenologically in the framework of Nearly Antiferromagnetic Fermi Liquid theory [|T3 
in some situations the short-range correlations are very strong and, therefore, a microscopic 
non-perturbative approach is needed to treat them. 

In the present Letter, we explore the Hamiltonian (|ID with half-fllled initial fermionic 
band by means of a trial wave function of Gutzwiller's type which, in addition to intrasite 
correlations, includes nearest-neighbor ones in an explicit form. Let us consider a system 
of N fermions of spin 1/2 on a lattice consisted of L sites. Then, in a general form, the 
correlated A^-particle trial wave function can be written as 



\^)=Il9x'\^o) (2) 

A 

where \ifo) is the A^-particle wave function of free fermions, for instance the Fermi sea 
IlkKkpt '^It ^kKkp, Q-fej. \0) , k IS the wave number of a fermion, kp^r is the Fermi wave number 
of fermions of spin a, g\ are the parameters which hes within interval ]0,cx3[, J\ are the 
projection operators onto all feasible configurations of a lattice site and a pair of adjacent 
sites. There are 4 such operators for itrasite configurations 
Xi = ^(l-nit)(l -nii),X2 = ^ni^{l-nii),X3 = ^ (1 - ni^) ni;,X4 = ^rai|ni|, 

i i i i 

(3) 
and 10 operators for nearest neighbor configurations, 

Yi = ^ (1 - riii) (1 - riii) (1 - nji) (1 - riji) , F2 = JZ ni^nnrij-^riji, and etc. (4) 

(see Table 1). 

The trial wave function (H) remains to be antisymmetric and conserves translational 
properties of Iv^o)- On the other hand, either intrasite and nearest-neighbor correlations are 
taken into consideration explicitly. From now on, we shall consider lattices for which the 
total number of nearest-neighbors pairs is equal to zL/2, where z is the number of nearest 
neighbors of a site. Let us denote normalized eigenvalues of the projection operators as 
xx 1$) = L^^Xx 1$), yx \^) = izL/2y^Yx |$). The eigenvalues turn out to be related to 
each other by normalization conditions 



5:x, = i, 5:/?Ai/A = i (5) 



and self-consistency conditions ||T^ 



yi + y3 + yi + y5 = xi, (6) 

1/2 + Z/3 + Z/8 + Z/9 = X4, 

2/4 + 1/6 + 1/7 + 1/8 = 3^2, 

1/5 + 1/7 + 1/9 + 1/10 = X3. 



As concentrations of fermions of each spins are fixed there are the only independent 
parameter xx and 7 independent parameters y\. In the case of half-filled band, additional 
constrains appear 

1/1 = 1/2, ye = yiO, 1/4 = 1/5 = 1/8 = 1/9 (7) 

which reduce the number of independent parameters yx to 3. It should be noted that, 
in contrast to (^ and (^, the condition (0) imposes the restriction on averaged values of 
yi only, that is, configurations, which violate (|^), enter into the trial wave function (^. 
Nevertheless, as will be shown below, one can omit them in the limit L — > oo. 

Assume that Xi = x^ = x, y^, y^, and y^ are independent parameters. Then, taking into 
account additional degeneracy of the half-filled band, the correlated trial wave function is 
transformed to 

1^) = gig^'gr'^'g^'^' l^o) (8) 

At first, we calculate the norm of the correlated trial wave function 

(^ I ^) = E W^.,ys,y.yr}gl'''gl'''''gT''V/''' = E R{^,y.y.,yr} (9) 

{x,yi,y4,y7} {x,y3,y4,y7} 

where the set {x, 1/3, 1/4, 2/7} describes configurations which contains Lx doubly occupied 
sites of the lattice, zLy^/2 nearest neighbor pairs of Y3 type and etc. A unessential factor is 
omitted in (|^) for simphcity. The summation is extended over all possible sets {x, 7/3, 7/4, y^}. 
^{x,ys,y4„y7} ^^ ^^^ number of configurations with the fixed set of the independent parameters. 
To calculate M/^{x,j/3, 1/4,5/7} we shall use the Kikuchi's pseudo-ensemble method |ll^ , |T5| . It 
should be mentioned that this method is exact on the Bethe lattice and approximated for 
lattices with closed paths [|14|. According to Kikuchi's hypothesis, we have 



where lower indexes oi W, T, and Q are omitted for the sake of simplicity. Q determines 
a number of arrangements of 10 elements corresponding to Yx taken zL/2 at a time. F 



is a fraction of correct arrangements in the pseudo-ensemble. In Eq.(^0]), dependent x\ 
and yx should be expressed in terms of x, 1/3, 1/4, and ^7. In the usual fashion, in the 
thermodynamic limit L -^ 00, we retain only the terms of the series @ which are very close 
to the largest one. As far as R{x,y3,y4„yr} is a positive function, one can search a maximum 
of its logarithm instead of the function itself. Let us retain the main terms on L only 
after taking the logarithm. It can be shown that this approach corresponds to substitution 
(zL/2)\ — > (L!)^'^ which was usually used |l^JT5[] . Then we obtain 



i^^2(.-i; 



2 J \2 



z [y2 In 2/2 + 2/3 In Vs + ^Va In 2/4 + 2/6 In y& + 2/7 In 2/7] 

(11) 



where 2/2 = a; — 2/3 — 2|/4, y^ = 1/2 — x — yj ~ 2|/4. The domain of function ([TT|) where we 
search the maximum is limited by conditions (^ and (|^). At its boundaries the gradient of 
function (^) is directed inwards the domain. That is why, the global maximum of In R/ L 
should be an inner one, and conditions din R/drjx = where rjx = x.y-^.y^^.yj are necessary 
for the global maximum. They lead to the following system of equations which relate Qi to 
X and yi 



X - 1/3 - 22/4 



(1/2 -x-yj- 22/4) (x - ?/3 - 2yi) 



2 2/7 

97- 



l/2-x-yT- 2y4 

To determine the expectation value of the Hamiltonian (j^) we need to calculate a density 
matrix of the first order: 

. {^\ E (ala,^ + H.c.)\^) 

"' = I WW) ■ ''" 

In contrast to GA, while a fermion hops from site i to j, it should be taken into account 
that configurations of adjacent pairs k — i and i — I change (FigJ^). Let us fix a particular 



lattice fragment (Fig.|T]), and calculate function W of the remain lattice from Eqs.(p^). 
Then, a fraction of configurations, which contain the fragment, can be found as 

TT// y{ij)Uyiki)Uyui) 

— = ^ '- . (14) 

Using (p!4D, we sum up all contributions to the density matrix. Then, it takes the following 
form 



Pi =4 



o / \2-l , Vsdr 2(2-1) , y79o93 2(z-l) 

2y4 {aia2) H a^ H 02 

9093 97 



(15) 



where Oi = (?/25'4 + y394/93 + 2/4 {97 + 1) /5'4) /x 

and a2 = {ye9i + y794/97 + 2/4 {93 + 1) /9i) /(1/2 — x). The first term of (|15D describes a 
motion of a fermion in the Hubbard subbands. The second and third ones arise from transi- 
tions between the subbands. In this expression, parameters go, g^, g^, g^ should be excluded 
by means of system (|12|). After some simplifications we find 



Pi = 8 (?/4 + y/ym) 



z-l 

' % + Vm + Vm + Vy7f 



_x{l/2-x) 
Finally, the total energy can be presented in the Gutzwiller's form 



(16) 



£=i.W^ = ,.„ + .!/ (17) 

where q = Pi/p?, p? is the density matrix (1T6|) at f/ = 0, ^o is the average energy of the free 
fermions. The ground state energy is determined as miuj^ ,^3 ^^ ,^^| (E). The function E turns 
out to be smooth, and its minimum is easily found numerically. It also should be noted 

that the total energy is invariant in respect of the following substitution: U — > —U,yi < > 

2/3)fl'7 ^ — ^ 93, 9o ^ — ^ l/S'o- A detailed discussion of the method presented above will be 
published elsewhere. 

As was mentioned above, the method used is exact for a Bethe lattice. We consider our 
results for z = 4 and 6 as approximate solutions for 2D and 3D lattices correspondingly. 
They are shown in the Fig.^, together with that of ID chain, and compared with the exact 
solutions for ID chain 0], D= 00 lattice [@], and numerical results of VMC method for 

6 



a paramagnetic phase |T^. The perturbative expansion (GA+1/D+l/D^) is not shown 



because it is very close to the GA curve for 2D and 3D (a shift of critical values of U 
between GA and GA+1/D+l/D^ is a few percent only ||12[)- To investigate the D= cxo 
limit of our model we carried out calculations for z = 50, 100, 200. They show that the 
ground state energy tends to the exact D= oo solution while z (or D) increasing. As one can 
see from the Fig.|^, at intermediate interaction of fermions ([/ ~ 1), our method gives the 
ground state energy significantly lower than that of VMC or GA+1/D+l/D^ procedures.. 
This means that nearest-neighbor correlation are essential in this region. 

Symmetric and antisymmetric correlation functions of the nearest neighbor 

Gs = {n^n^y + {n^n^y = 2 (1/2 + 2i/4 + Ve) , (18) 

Ga = (^T^i)' + (^i^t)' = 2 (1/2 + 2|/4 + y?) 

are shown in Fig.^ for the same lattices as in Fig.|^. The prime in Eqs.|l^ denotes the 
averaging over nearest-neighbor pairs only. One can see in Fig.^ that even at f/ = some 
correlations appear due to a exchange hole. An increasing of U leads to enhancing of 
correlations but there exists a saturation point. At further increasing of U the nearest- 
neighbor correlations remain almost constant and variation of the ground state energy is 
due to intrasite correlations only. 

In conclusion, a new non-perturbative approach to problem of strongly correlated 
fermions is reported. A trial wave function which includes nearest-neighbor correlations 
is constructed. For a half-filled initial fermionic band, the ground state energy of the Hub- 
bard Hamiltonian as well as correlation functions are calculated for 2D and 3D lattices and 
the ID chain. Since the correlated wave function \ip) has the same translational properties 
as \(po) does, all the results obtained above describe a paramagnetic state. To build entire 
phase diagram for Hubbard model our approach need to be extended to an antiferromag- 
netic phase. The method may be especially beneficial when systems with strong short-range 
correlations, for example CUO2 planes of HSTC, is considered. 
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FIG. 1. A fragment of a tree {z = 4). While a fermion hops from j to k site, configurations of 
adjacent pairs need to be taken into account. 

FIG. 2. The ground state energy obtain by minimization of Eq.(17) for ID chain (z = 2), z = 4 
(2D), and z = 6 (3D) lattices (solid lines); the exact solutions: GA ( dotted line) [2], ID chain 
(dash-dot line) [4]; and numerical VMC calculations (dash line). 

FIG. 3. The correlation functions Gg (dashed lines) and Ga (solid lines) of ID chain, z = 4 (2D) 
and z = 6 (3D) lattices. 

TABLE I. Pair projection operators, corresponding configurations and the degeneracy factor 
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